The xmsinterp library contains classes for performing 2d (x,y) linear, natural neighbor, and idw interpolation: InterpLinear and InterpIdw. Natural neighbor interpolation is performed using the linear interpolation class. 3d (x,y,z) interpolation can be performed using the idw class. Points (x,y) must be given to create the interpolation class.
We will learn how to use the interpolation classes by using a set of 2D points (x,y) with a measurement at each point. These points come from a csv file. Then we will generate a grid that bounds these points and use different interpolation options to visualize the data.
import param, panel as pn, holoviews as hv
from holoviews.streams import Pipe
from bokeh.document import Document
import numpy as np
import xmsinterp as interp
import pandas as pd
import xarray as xr
%matplotlib inline
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import time
hv.extension('bokeh', 'matplotlib')
# make the plots bigger by default
plt.rcParams['figure.dpi'] = 120
# read in the data to do the interpolation
df = pd.read_csv("plumedat.csv")
plot = df.plot(kind='scatter', x='x', y='y', c='c', title='sampled concentration', cmap='rainbow')
# drop the id columns from the data frame so we can pass the remaining data as points to the interp class
df.drop(['id'], axis=1, inplace=True)
Get the min, max (x, y) of the data and generate a grid of points.
# get the (x,y) range of the data and expand by 10%
x_range = (df.loc[df['x'].idxmin()]['x'], df.loc[df['x'].idxmax()]['x'])
buf = (x_range[1] - x_range[0]) * .1
x_range = (x_range[0] - buf, x_range[1] + buf)
y_range = (df.loc[df['y'].idxmin()]['y'], df.loc[df['y'].idxmax()]['y'])
buf = (y_range[1] - y_range[0]) * .1
y_range = (y_range[0] - buf, y_range[1] + buf)
# adjust these values for a higher/lower resolution grid
number_of_x_points=600
number_of_y_points=400
x = np.arange(x_range[0], x_range[1]+1, (x_range[1]-x_range[0])/number_of_x_points)
y = np.arange(y_range[0], y_range[1]+1, (y_range[1]-y_range[0])/number_of_y_points)
grid = np.meshgrid(x, y)
pts = [x for x in zip(*(x.flat for x in grid))]
xx, yy = grid
Interpolate to the grid using linear interpolation and display the results with matplotlib.
linear = interp.interpolate.InterpLinear(df.values)
z = linear.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
# make the plots bigger by default
plt.rcParams['figure.dpi'] = 120
p = plt.contourf(xx,yy,zz)
cbar = plt.colorbar(p)
Here is an example visualizing with xarray.
xa = xr.DataArray(np.reshape(linear.interp_to_pts(pts), (len(y), len(x))))
xa.plot()
Interpolate using natural neighbor. The first example shows a "constant" nodal function and the second example shows a "quadratic" nodals function. A nodal function can better model trends that exist in the measured values.
linear.set_use_natural_neighbor()
z = linear.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz)
cbar = plt.colorbar(p)
Now do natural neighbor with a quadratic nodal function.
linear.set_use_natural_neighbor(nodal_func_type="quadratic", nd_func_pt_search_opt="nearest_pts")
z = linear.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz)
cbar = plt.colorbar(p)
Now we will use Inverse Distance Weighted (IDW) to interpolate to the grid. Notice that unlike Linear and Natural Neighbor interpolation, IDW interpolation can extrapolate beyond the bounds of the input data points.
idw = interp.interpolate.InterpIdw(df.values)
z = idw.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz)
cbar = plt.colorbar(p)
Now interpolate using a nodal function to try to better capture trends in the data set.
idw.set_nodal_function(nodal_func_type="gradient_plane")
z = idw.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz, 10)
cbar = plt.colorbar(p)
Notice, in the previous interpolation that we get negative values in some locations. If we assume that there should be no negatives values (like when measuring concentration) we can use the truncation options to adjust the interpolation.
idw.set_truncation_max_min(smax=150, smin=0)
z = idw.interp_to_pts(pts)
zz = np.reshape(z, (len(y), len(x)))
p = plt.contourf(xx,yy,zz, 10)
cbar = plt.colorbar(p)
Below we define a param class for the interpolation options. Then we use panel to show widgets and a dynamic plot. The user can change the widget values and then push the "Compute new plot" button to update the plot using the currently specified interpolation options.
pipe = Pipe([])
bounds = x_range[0], y_range[0], x_range[1], y_range[1]
class InterpOptions(param.Parameterized):
interpolation_option = param.ObjectSelector(default='idw', objects=['linear', 'natural_neighbor', 'idw'],
precedence=1)
nodal_function = param.ObjectSelector(default='constant', objects=['constant', 'gradient_plane', 'quadratic'],
precedence=1.1)
truncation = param.Boolean(default=True, precedence=2)
truncation_maximum = param.Number(default=df.loc[df['c'].idxmax()]['c'], precedence=2.1)
truncation_minimum = param.Number(default=df.loc[df['c'].idxmin()]['c'], precedence=2.2)
compute_new_plot = param.Action(lambda x: x.update_hv_plot(), precedence=0)
index = 0
@param.depends('interpolation_option', watch=True)
def _update_nodal_function(self):
if self.interpolation_option == 'linear':
self.params('nodal_function').precedence = -1
else:
self.params('nodal_function').precedence = 1.1
@param.depends('truncation', watch=True)
def _update_truncation_opts(self):
if self.truncation:
self.params('truncation_maximum').precedence = 2.1
self.params('truncation_minimum').precedence = 2.2
else:
self.params('truncation_maximum').precedence = -1
self.params('truncation_minimum').precedence = -1
def generate_zz(self):
# Get Interp Object
interp_object = None
if self.interpolation_option != 'idw':
#interp_object = linear
interp_object = interp.interpolate.InterpLinear(pts=df.values)
if self.interpolation_option == 'natural_neighbor':
self.index += 1
interp_object.interp_to_pt((0,0)) # this will force triangle creation
interp_object.set_use_natural_neighbor(nodal_func_type=self.nodal_function,
nd_func_pt_search_opt="nearest_pts")
else:
interp_object = interp.interpolate.InterpIdw(pts=df.values)
interp_object.set_nodal_function(nodal_func_type=self.nodal_function)
if self.truncation:
interp_object.set_truncation_max_min(self.truncation_maximum, self.truncation_minimum)
z = interp_object.interp_to_pts(pts)
bob = np.reshape(z, (len(y), len(x)))
return bob
def update_hv_plot(self):
pipe.send([self.generate_zz()])
def get_image(self):
return hv.Image(self.generate_zz(), bounds=bounds)
opts = InterpOptions(name="")
#opts.get_image()
%%output size=200
dmap = hv.DynamicMap(lambda **kw: opts.get_image(), streams=[pipe])
pn.Row(opts, dmap)